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Abstract 

The phase diagram at zero temperature of a lattice SU(2) scalar-fermion 
model in (2+1) dimensions is studied numerically and with Mean-Field meth- 
ods. Special attention is devoted to the strong coupling regime. We have 
developed a new method to adaptate the Hybrid Monte Carlo algorithm to 
the 0(3) non-linear a model constraint. The charged excitations in the vari- 
ous phases are studied at the Mean-Field level. Bound states of two charged 
fermions are found in a strongly coupled paramagnetic phase. On the other 
hand, in the strongly coupled antiferromagnetic phase fermionic excitations 
around momenta (±^,±5-,±^) emerge. 
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1 Introduction 



The model we are going to study was proposed in references Q and |||] as a natural 
extension of the lattice 0(3) non-linear cr-model in (2+1) dimensions to include charge 
carriers. It is a lattice model of interacting spins and Dirac fermions in (2+1) dimensions, 
with only two free parameters in addition to the temperature: a nearest-neighbour spin 
coupling and a spin-fermion coupling. The model describes quantitatively some of the 
features of the doped copper oxide compounds |l|, ||. 

In the present article we want to present a careful, detailed discussion of the model, its 
symmetries, and its properties, and give full technical details and results of the Mean-Field 
(MF) and Monte Carlo (MC) calculations, some of which were reported in Ref. jy]. In 
this paper, our Mean-Field and numerical studies will be limited to the zero-temperature 
case, corresponding to infinite Euclidean time direction. 

The remainder is laid out as follows. In Section [2] we present our model, discuss 
the choice of lattice fermions, comment on the symmetries of the model, give its phase 
diagram and prove the reality of the fermion determinant, even in the presence of a 
chemical potential. In Section ^ we examine the phase diagram of the model in the MF 
approximation. Our MF calculations are based on small- and large-?/ expansions combined 
with saddle points methods. The method allows us to handle (products of) fermionic 
variables occurring in the expansion of the fermion determinant in a well-defined way. 
In Section |I| we use MC simulations to complete the study of the phase diagram. For 
this purpose we have developed a new method that exactly solves the technical problem 
related to the length-1 constraint on the spin variables ||. Section || is devoted to a 
study of the relevant excitations in the different phases of the system, at the MF level. A 
crucial result is the dynamical generation of spin singlet bosonic bound states of charged 
fermions in the so-called paramagnetic strong (PMS) phase. At the MF level we have not 
found evidence for fermionic excitations at zero temperature in this PMS phase. Another 
interesting result is the emergence of fermionic excitations around momenta (±f , =t§, ±f ) 
in the strongly coupled Antiferromagnetic ( AFM) phase 0] . Finally, Section ^ is devoted 
to our conclusions and projects. 



2 The model: Formulation, Symmetries, Phase 
Diagram 



The model is defined by the following (2+l)-dimensional lattice Euclidean (imaginary 
time) path integral, 



We use this expression as our starting point, but it should be noted that the model depends 
only on the ratio y = X/p, through a change in the normalization of the fermion field. In 




(1) 



with action 



S = ~ (j) x ■ <f> x+[l + -j^xl^i^x+ji ~ if> x -fi) + Yl ■ • (2) 
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terms of the effective spin-fermion coupling y, we get: 



s=-^2k(j> x - (f> x+fl + ^2^ v^Ox+a - $r-/0 + Yl y ^ x ^x • • (3) 

X,fJ, X,fJ, X 

Here x runs over a (2 + l)-dimensional cubic Euclidean space-time lattice. Each tp x is 
a fermionic four-spinor as a shorthand for two flavours of two-component Dirac spinors. 
Both flavours are taken in the same irreducible spinor representation, with 2x2 gamma 
matrices taken as the Pauli matrices a^. The 4x4 matrices 7^ in Eq. (||) have the form 

T" = {H ^ =1 ' 2 ' 3 - (4) 

The kinetic term for the fermions is of the nearest-neighbour (hopping) form. Lattice 
fermions defined in this way undergo species doubling in the perturbative continuum 
limit ||. For two reasons we are going to leave this matter aside in this work. First, 
we are particularly interested in the strong coupling non-perturbative regime where more 
of the interesting physics occurs In this strong coupling regime all the fermions, 

the physical one as well as the doublers, decouple in the continuum limit ||. Second, this 
model described qualitatively some of the features of the doped copper oxide compounds 
||, where the lattice space is given by the material. 

The three-component are real scalar bosonic variables, subject to the constraint 
4> 2 = 1, as in the 0(3) non-linear a- model. The last term in Eq. @ describes the 
interaction between cf> and the Dirac fermions, which is diagonal in Dirac space. The 
Pauli matrices r a act in flavour space. 

Let us now consider the symmetries of (0). First of all, we have the usual cubic 
symmetry. Next, there is a discrete parity symmetry, which in (2+1) dimensions is defined 
as the reflection of one of the spatial axes, say the x-axis. Under this parity symmetry, 
the fermions can be seen to transform as 

ip^a-Lip, tp^-tpai, (5) 

so (f) is a pseudoscalar in this sense. In addition, the action (|3|) is invariant under an 
SU(2) flavour symmetry in which ip transforms as the fundamental representation and 
transforms as the adjoint one. Note that by requiring the two fermion flavours to have 
the same Lorentz structure (that is, by choosing the 7's as in @) no fermion mass term 
is allowed which preserves the above symmetries H. 

There are two more discrete symmetries of our model (||), which will be useful in the 
MF calculation of the phase-diagram. The first one is 

Z(k,y) = Z(k,-y), (6) 
which becomes clear if we make the change of variables 

tp x -> exp (i^xj ipx, ipx ->• exp y^xj (7) 

where 

= (8) 
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This implies that Z(k, y) is a function of y 2 only and we can restrict ourselves to y > 0. 
In addition, there is a symmetry 

Z(k,y) = Z(-k,-iy), (9) 

as can be seen by making the substitutions 

ip x -> exp (ij^xj ipx, ipx ->• exp V^, 0z e x . c/v (10) 

The latter symmetry implies that the lattice regularization of the non-linear cr-model, 
y=0 (or y=oo, see Sections ||, ||), is equally valid in a ferromagnetic or an antiferromagnetic 
phase. 

In order to perform computations in models of this type, one has to integrate out the 
fermions. This integration leads to a 0-dependent fermion determinant. It is important 
to know whether this determinant is a real number. To study this, let us write down 
the original fermion matrix (Latin letters x,y, . . . will refer to lattice points, i,j,..., will 
represent flavour indices, while Greek letters a, /?,... are used for Dirac indices): 



M X ai;yf3j — K-xai;y0j ~r" Y xa i-yf3j j (H) 
K X ai\yf3j = ~ ^ / (^+A t .y ~ ^x—fj,,y) (12) 

Y xa i;y/3j = y$xy{4>- T )ij (13) 

Keeping in mind that for Pauli matrices 020^2 = — cr*, where * means complex conjuga- 
tion, and that [7, r] = 0, one easily proves that, for real y, 

a 2 r 2 (K + Y) g 2 t 2 = — (K* + Y*) . (14) 

Therefore, 

det (K + Y)= det (-K* - Y*) = [det (K + Y)]* , (15) 

i.e. the determinant is real. Thus, by doubling the number of fermion families, we obtain 
a positive determinant. Had we introduced a chemical potential, /j,, the only change would 
be the introduction of e^ on the temporal links of the kinetic matrix ||. The essential 
requirement for Eq. (O) to hold (that the only non-real numbers are in 7,r) is thus not 
endangered by the chemical potential and the determinant is still real. 

The phase diagram of the model at zero temperature is shown in fig. |]. Notice that it 
is very similar to the phase diagram of (chiral) Yukawa models for the electroweak sector 
of the Standard Model of elementary particle interactions Q. At y = 00 and at y = we 
recover the non-linear cr-model (see sections ^J2|) with its well known paramagnetic (PM), 
ferromagnetic (FM) and antiferromagnetic (AFM) phases. At finite y, we expect these 
phases to extend into the (k,y) plane. One of its most important features is that there are 
two mutually disconnected paramagnetic phases, one at weak coupling (called PMW) and 
one at strong coupling (PMS). One sees that the PMW-FM and the PMW- AFM transition 
lines meet in a point A, where this disordered phase ends. In the strong coupling sector of 
the phase diagram, a similar behaviour is found, with the two transition lines meeting at 
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Figure 1: Phase diagram of the action (||), for two fermion families. Dashed lines 
are from the MF calculation, solid lines from a MC calculation on an 8 3 lattice. 

point B. This observation means that one may expect totally different behaviour in each 
of the two paramagnetic phases. This is indeed the case, as we shall see later. 

As there is no evidence for a phase transition between the strong- and weak-coupling 
regions of the FM and AFM phases, we name them FM(W) and FM(S), AFM(W) and 
AFM(S) (note the parentheses). There may be crossovers between these regions, though. 

Between the points A and B, we find a phase where both the magnetization and the 
staggered magnetization are different from zero. We name this phase exotic magnetic 
(EM). An appealing possibility is that it corresponds to a helicoidal phase. We expect the 
EM phase to disappear for large enough —k, but we have not explored this numerically. 



3 Mean Field Calculations of the Phase Diagram 

Our aim in this section is to determine the zero-temperature phase diagram of the model 
in the y-k plane (cf. Fig. [|), using Mean-Field techniques. These calculations already 
provide a lot of insight, especially for the strong coupling region. They will be contrasted 
with numerical simulations for the phase diagram in Sect. ||, and they will be extended to 
a study of the relevant charged (quasi-particle) excitations in Sect. [|. 

Our Mean-Field calculations are based on small- and large-y expansions combined with 
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the saddle point methods described in Ref. [1C]. This approach guarantees a systematic 
expansion in 1/d, which is particularly important for operators which are zero to lowest- 
order. Our particular method furthermore allows us to handle (products of) fermionic 
variables occurring in the expansion of the fermion determinant in a well-defined way. 
These techniques have been developed and applied in the context of similar lattice models 

12] of the Electroweak sector of the Standard Model of elementary particle interactions, 
and in the study of the antiferromagnetic (ft 4 model ]l3|] . 

We shall first concentrate on the small-y region, and incorporate the fermion determi- 
nant up to 0(y 2 ). 

In order to apply the saddle-point method, the integration over the fields must be unre- 
stricted. We therefore need to replace the integration over the spin vectors eft, constrained 
by the condition \cj)\ = 1, with an integration over unconstrained variables This is done 
by multiplying the functional integrand in Eq. (jl|) by 

3 
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and replacing a conveniently chosen subset of the eft variables in the action S with £ fields. 
We obtain 



Z 



D£DA 
T ex P 



^ ' L 2;, A* x 

DtftDtft exp - ^2 ^x^i^x+fi ~ ^x-fx) 

n | / exp ^ iAx ' ^ - v i'x 4> x ■ Tipx] 



(16) 



Note that both the £ fields and the auxiliary fields A are unconstrained. 

Now we have to integrate out the constrained variables (ft^ (as well as the fermions), 
before the mean fields can be introduced. Let us concentrate on a single cft n integration 
in Eq. (|l6|), dropping the subscripts n for simplicity. First, we perform an expansion in 
powers of y. We can write 



— exp [iA ■ eft - y $ (ft ■ Tift] 



= exp [u(iA)] exp 

where we have defined 
Q a = {ft T a ift , u (iA 

and we have introduced the notation 



yQ a -m iA + ^y Q a Q T +°(v 



(17) 



In /^exp[iA-0] , T ab 



'iA> 



(0) iA = y ^ O exp [iA ■ (ft] j |g exp [iA 
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In addition, we introduce a Hubbard-Stratonovich vector parameter A to deal with the 
quartic fermion term in Eq. (|i~7|), 



cxp 



\y 2 £ Q a Q b T ab 



a.b 



dX 



cxp 



x exp 



(2tt) 3 /2 



- y A a A a 

2 ^ 

a 

g a A 6 



(18) 



(Note that the matrix T is self adjoint, and positive definite if A is imaginary, so the 
square root is well defined). Thus, up to this order in y 2 , the action is bilinear in the 
fermion fields. 

Carrying out the fermion integration in Eq. (|l6|) now gives det M, where 



M 



x,a,i;y,f3,j 



K x , a , iwAj + yS xy S aP ^((Oa, - £ (V^T A 6 )r§ ■ 



(19) 



The matrix K has been defined in Eq. (|12[). 

The mean fields are the field values at the saddle point of the free energy 

-F = j2HiA x ) + kj2z x -Zx+„-iJ2 A *-z*-lJ2 x l +TTlo z M - 



(20) 



x,ii 



A choice of the mean fields should be done at this point, as we cannot calculate log Af for 
general {A x , X x }. An appropriate choice for the study of a PM-FM phase transition is 



X x 

in terms of which (N is the lattice volume) 



(0,0, -ia) 

(o,o,«), 

(0,0, A), 



(21) 



F/N 



-u[a) 



kdv 2 + av H — A 2 Tr 



N 



logM, 



with a, u and A satisfying the saddle point equations 



VF 



(a,v,X) 



0. 



(22) 



(23) 



The fermion matrix, M(a,v, A), can be calculated in momentum space, where it is 
diagonal in its momentum indices. One easily finds 



det M = exp 



2j> 



Eu=ism P M 



(24) 
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where we have divided out the determinant for free fermions. We need only the leading 
0(y 2 ) contribution to the exponent, hence the mean field free energy becomes, in the 
infinite volume limit: 

F/N = -u(a) - kdv 2 + av + ^A 2 - 2y 2 (u'(a) - Xy/u"(a)^j 2 C , (25) 

where 

f 71 d 3 v 1 

C = / 777-^^3 —~ — = 1.0109240. (26) 



U(2vr)3 Y^U^p, 

Next, we shall discuss the actual solutions to Eqs. (f23|). From u(a) = In (sinh a/a), 
one easily finds that a = v = A = always fulfill them. For small k, y, it is a true 
minimum of the free energy. This characterizes a paramagnetic (PM) phase, since none 
of the fields develops an expectation value. 

For larger values of k and y, there is another, non-trivial solution, corresponding to a 
ferromagnetic (FM) phase. It emerges when a negative mode in F/N starts to develop, 
as a function of the mean fields, and the transition between the two regions is given by 
the condition (F" is the Hessian matrix) 

detF"L _ n , n , = 0. (27) 
This condition is satisfied for F/N of Eq. (p5|) if 



3 2 C() r> 

"s—r*- (28) 

This curve defines the phase transition line between the PM and FM phases in the 
small- y region. Using the symmetry (0), we deduce that there is a similar transition 
separating the PM and AM phases, 

Let us finally remark that in the presence of Nj such fermion fields we would have Nf 
factors of det M, leading to a multiplication of Co by Nf in Eqs. (pSQ and (p9|). 

The large-y region is easier to deal with. Here it is convenient to integrate out the 
fermions directly in Eq. (|l6|), leading to (summation over repeated index is carried) 

det M Xj0 , ti .y^j = det ^K x>aji . y> pj + yb afi b xy ^ <f x rf^j (30) 
= det (^yS ai 5 xz ^2^T^j 

X det ^SzyS-ypSkj + - 4 >b x T kl K z,-(,l;y,f3,j S j • 



(31) 



Here we have used that (^ a <ft a " ra ) 2 = 1 (recall the 0's are unit vectors). Now we can 
expand log(det M) in powers of 1/y. The 0(l/y) term vanishes by virtue of K xx = 0. To 
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second order one obtains 

log det M = log y AN + Tr £ <f> a x r^K xai]tll £ ^K^^j (32) 

= \ogy AN + \ ^cj) x -(f> x+fM . (33) 

Here, log y AN is an irrelevant constant that can be dropped. Notice also that this expression 
will acquire a prefactor Nf if there are Nf identical fermion flavours. One sees that, up 
to 0(l/y 2 ), the only effect of the fermion determinant is a renormalization of the scalar 
hopping parameter of the 0(3) model, 

k -» k + N f \ . (34) 

Note that we did not introduce any mean fields to derive this result. The usual MF 
treatment of the 0(3) model with this renormalized coupling now immediately gives us 
the required phase transition lines in the large-y region of our model: 



±o3 - N f -2 ■ ( 35 ) 



"2d V 

It is interesting to compare the small- and large- y results, to leading order in 1/d. As 
is well known, the first order in this expansion is equivalent to any MF approximation, up 
to higher-order terms. For this purpose, we need the 1/d expansion of the constant Co in 
Eq. (^), which can be calculated as follows: 



oo 



Co{d) = Fi^ ^ l — =2 I ds ^~ SI ^ d ( 36 ) 

-3K^(?))- 

where Io(s) = (d8/2-K) exp(scos#) is the modified Bessel function. In fact, the second 
equality in Eq. ( |36l) was used to obtain the numerical result (|26|) for Co- 

Keeping only the leading-order term 2/d for Co we find that the phase transition lines 
would meet at y 2 = 2/d. 

Now we are ready to map out the phase diagram of the model, as predicted by the MF 
method for the weak and strong coupling regions. This is done in Fig. [l]. The vertical axes 
at y = and y = oo correspond to the 0(3) model, with its disordered (PM) and ordered 
(FM and AFM) phases. These phases extend into the y-direction, both for y > and 
y < oo. Note that all the phase transition lines bend downward. This can be understood 
intuitively by assuming a MF value for the fermion condensate, which would act as an 
external field tending to align the spins <f> in parallel. 

4 Monte Carlo: Method and Results 

A well established method for dynamical fermion simulations is Hybrid Monte Carlo 
(HMC) jl4|]. However, the implementation of this algorithm in a model with constrained 
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variables is not straightforward. This has been satisfactorily achieved for models with vari- 
ables belonging to a Lie group [15], like SU(iV) gauge theories or like some spin-models, 
such as the 0(N = 2,4) non-linear a-models. However, for other spin variables (not in a 
Lie group), as in the 0(3) non-linear cr-model, this had not been satisfactorily solved yet, 
although the problem arose already in the first simulations using the Langevin algorithm 
||. Our solution is a generalization of the strategy in 

We shall first discuss our solution in the quenched approximation, where comparison 
with other algorithms is possible (Section 4.1), and then deal with the full theory in 
Section [4.2[ Finally our Monte Carlo results for the phase diagram of the full theory will 



be presented in Section 4.3. 



4.1 The HMC method for the quenched approximation 

For the purpose of discussion it will prove convenient to briefly describe the HMC method 
for unconstrained bosonic variables 4>(x) with action Sb^) (see ref. |]l6| for a pedagogical 
presentation): 

1. Introduce uncorrelated gaussian variables ir(x) of unit variance (the conjugate mo- 
menta for the fields (f>) and define a Hamiltonian 

H = Y j \^)+Sb{^). (38) 

X 

One can then use the hamiltonian equations of motion 

tt(x,t), (39) 
_SSb_ 
54>(x,t) ' 

to perform a microcanonical Molecular Dynamics evolution in "Monte Carlo time" , 
t. After a certain period of MC time (called "trajectory"), new random momenta 
tt(x) are chosen ("refreshing" the momenta). The crucial properties of Eqs. ( f40| ) are 
their time reversibility, and the invariance of the Liouville measure, Deft Dir, under 
time evolution. 

2. In practice, the molecular dynamics equations of motion for a trajectory are dis- 
cretized into N steps At. This is done using a leap-frog algorithm which is exactly 
time reversible, but does introduce a systematic error which shows up as a non-zero 
AH = 0(At 2 ). The detailed-balance is not endangered by this error, because a 
Metropolis acceptance step is performed. For fixed trajectory length, N can then 
be tuned to optimize the overall efficiency. 

To generalize the method to constrained variables, one needs to appropriately define 
the conjugate momenta and the equations of motion in order to preserve the constraint 
and, most importantly, not to spoil the time reversibility. Each spin variable, <j>, lives on 
the surface of a two-sphere, and correspondingly one could imagine an algorithm with two 
independent conjugate momenta, maybe living in the perpendicular plane (</> • 7r = 0). 
However, changing the constraint from the field <f> to the momenta is not very appealing 



tt(x,t) = 
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(and, from the practical side, one would need to worry about two constraints in the numer- 
ical integration). A different approach, the use of spherical coordinates, has the drawback 
of a non-planar integration measure. Our very simple algorithm avoids constraints and 
non-planar measures, by introducing three conjugate momenta per spin. 

We shall start from an analogy with the dynamics of a particle living in the sphere, a 
potential (V") acting on it. The Hamiltonian is 

T 2 

^sphere = _ + y^y ^ 

Here L is the angular momentum, 0x0. The equations of motion can now be obtained 
from the Poisson Bracket with the hamiltonian (flC|): 

= Lx0,L = -0x^. (41) 

In this expression J3£ stands for (jfa, 

This formalism is still inconvenient for us, because the constraint • L = complicates 
the generation of random momenta according to a Gaussian distribution. However, the 
following simple facts can be straightforwardly established from the equations (f4l|): 

I. Both 2 and • L are conserved through the time evolution. If the initial condition 
verifies the constraints 2 = 1 , • L = 0, this will not be spoiled by the dynamics. 

II. The dynamics is time-reversible. 



III. Although the Li cannot be all canonical variables 17], the "Liouville" measure, 
D(f) DL{= d(f>i dcp2 d(p3 dL\ dL2 dL%), is left invariant by the time-evolution. 

IV. The Hamiltonian is a constant of the motion. 

Now let us forget about the constraint 0- L = 0, i.e. we introduce a new field P which can 
have a "radial component" (it is no longer an angular momentum), but we keep the Eqs. 
of motion (|4l|). Obviously, statements I-IV will still hold. Whether a symplectic structure 
is hidden under this new dynamical system is unclear, but also irrelevant (properties II 
and III are the essential ones for HMC to be a correct algorithm p6[ ). 

So, we introduce three momenta per spin, P = (Pi, P2, P3), and write down the 
Hamiltonian 

H = ^^! + 5 B (0). (42) 

X 

Equations of motion respecting properties I-IV are easily generalized: 

S S b 

4>(x,t) = P (x,t) X 4>(x,t) . P(x,t) = -4>(x,t) x • ( 43 ) 

As expected, the evolution equations for the S 2 fields cj) take the form of (infinitesimal) 
rotations, while the conjugate momenta can be considered as living in the Lie Algebra of 
SO (3). The discretized leap-frog form of these equations is therefore naturally formulated 
in terms of finite SO(3) rotations, 

x (nAr + Ar) = exp[ATP x ((n + -)Ar) ■ J] K (nAr) , (44) 

P x ((n + hAr) = P x ((n - I) At) - (a; , nAT) x J Sb At , (45) 
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Table 1: Values for several observables in the quenched model (|3]) on an 8 3 lattice at 
k = 0.693 ~ k c , obtained with our implementation of HMC and with Wolff's single 
cluster algorithm 



Algorithm 


(E) 


d k (E) 


x/v 


£ B 


HMC 


0.3505(5) 


1.51(2) 


0.1426(9) 


4.47(2) 0.800(6) 


Wolff 


0.35061(13) 


1.501(10) 


0.1432(2) 


4.486(9) 0.8031(18) 



where J are the 3x3 generators of SO(3), satisfying 

(exp[6>n • J])ij = 5ij cos 8 + riiUj (1 — cos#) — e^fc^-fc sin 9 (46) 

for unit vectors n. Again, the length constraint on the eft fields is preserved by construction. 

This final result is reminiscent of the elegant solution for models with variables be- 
longing to a Lie group and conjugate momenta in the group algebra (or vice versa) fl5|| . 

In our case, «Sb quenched = — k ^ n 4> n ■ </> n+ ^, so the HMC algorithm can now be 
implemented in a straightforward manner. To test the algorithm, we have simulated 
the 0(3) model on an 8 3 lattice at k = 0.693 ^ k c @ with our HMC algorithm and 



with Wolff's single-cluster embedding algorithm [18]. Let us first define the measured 
observables, and then compare them. 

In this work we have only measured bosonic observables, as our sole objective was the 
numerical determination of the phase diagram. We have constructed our observables in 
terms of the Fourier transform of the spin field: 

X 

where V = L 3 is the lattice volume. 

We define the non-connected finite-volume susceptibilities as 

X = ^(m 2 (0,0,0)>, Xs = y<m 2 (7r,vr,7r)>. (48) 

The subscript 's' on Xs stands for 'staggered', and this term is used to label quantities 
which are taken with a weight —1 for the odd lattice sites, corresponding to momentum 
(tt, 7t 3 7r). Notice that xjV is a pseudo order parameter, which should be of order one in a 
ferromagnetically broken phase, and of order 1/V in a paramagnetic or antiferromagnetic 
phase (and similarly for Xs/V)- 

Another quantity of interest is the Binder cumulant 

5 3 ((™ 2 mo)y 

2 2 (m 2 (0,0,0))^ 

with an analogous definition for the staggered variant B s . 

One expects B = 1 in the FM phase, where x/V is non-vanishing in the thermodynamic 
limit, while it should be of order 1/V in the PM phase, far from the phase transition. 



B = tt- ^ ^ , (49) 
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For the correlation length, we use a definition which is easy to measure and gives 
accurate results: 

^{Sm)' 2 ' (50) 

where F is the squared Fourier transform at minimal non-zero momentum, 

F = — ^|m(2vr/L,0,0)| 2 ^> + permutations) . (51) 

Again, the generalization to staggered quantities is straightforward. Another kind of 
observable, needed for the standard extrapolation method [^] , is the normalized nearest- 
neighbour energy 

E = v^-^+A> = ^Z. (52) 

We also measure its fluctuation, given by 

3V (<£ 2 > - (Ef) = ± (E) . (53) 

In Table |l] we compare the values obtained for these observables, using our HMC 
algorithm and the single-cluster algorithm. We find excellent agreement. Of course the 
efficiency of our implementation of HMC is not competitive with a cluster method in the 
0(3) non-linear a- model. But it could be useful in other models where cluster methods 
are not effective in reducing the dynamical critical exponent z (for instance, when some 
kind of frustration is present fl2l"| ), while HMC is expected to yield z = 1 for any bosonic 
model. 



4.2 The HMC algorithm for the full theory 

The only restriction imposed on HMC is that the fermion bilinear in the action should be 
given in terms of a positive definite matrix. This will the case if we consider two identical 
fermion families (Nf = 2) as is usually done in lattice gauge theories. After integrating 
them out we obtain (detM) 2 = det (M^M), where M is the fermion matrix for a single 
fermion family. As we are mainly interested in the strong spin-fermion coupling region, it 
makes sense to perform the following manipulation: 

det M = det {Y + K) = y 4V det (1 + Y~ X K) (54) 

(cf. Eqs. (30,31])). The constant factor y iV can be dropped, and we define M = l + Y^K. 

Next, one re-exponentiates the (inverse) fermion matrix by introducing the so-called 
pseudo-fermions z x , which are complex four-component c-number fields. The partition 
function is then 

Z = J Dcf)DzDz exp (^—Sb — z(M^ M)^ 1 z^j . (55) 



For further details we refer to Ref. p6 |. 



Now the HMC Hamiltonian becomes 



H = E \ p l - k E ^ • ^ + z] ( MtM ) 1 z > ( 56 ) 



X,fl 
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and the time reversible, constraint and energy preserving equations of motion are 



0(as,T) = P (x,r) X 0(as,T)» ( 57 ) 
■P(x,t) = ~ k Yl {^(x+^t) + 0(x-m,t)) X 0(*,t) 



x 4>{x,r) m + h - c - 

x,t) J 



-i 



z. 



For the inversion of the fermionic matrix, we have employed the conjugate gradient 
algorithm. To formulate the stopping criterium, let us define h = [M*M s j z, h n being 
the n th trial solution. We continued the conjugate gradient iteration until 



h I 2 
'ml 



In the simulation, we need the inverse matrix both for the leap-frog and for the 
Metropolis accept-reject step. It is clear that R does not need to be the same in both 
cases. For the Metropolis step, lack of accuracy in the inversion will bias the simula- 
tion. To control this, we have checked that the Creutz parameter (exp(— AH)) equals 1 
within errors. In some regions of parameter space R values as small as 10 -25 were needed. 
The essential requirement on the leap-frog is full reversibility in the numerical integration 
of the equations of motion (up to the numerical precision reachable with 64-bit floating 



point arithmetic). As first noticed in ref. 22], this has no relation with R if the seed 



for the conjugate- gradient algorithm is chosen to depend on the actual configuration only 
(ho = z, for instance). However, if R is too large, the numerical integration will produce 
large changes in the Hamiltonian, and the Metropolis acceptance will be poor. We have 
found that R = 10~ 7 during the leap-frog steps allows for a 50% acceptance. 

In a first implementation of a new MC algorithm, some consistency checks are ex- 
tremely useful. In addition, there are three parameters to be adjusted for optimal perfor- 
mance, N,t and R. We have carried out the following tests: 

1. We have explicitly checked reversibility of the leap-frog algorithm. 

2. We have checked that (exp(-AH)) = 1 within errors. 

3. The gaussian expectation values, /z^ (M^M) 1 z^ = 4 and (-P 2 ) = 3 have been 
checked. 

4. We have checked that AH oc (At) 2 in the leap-frog integration, for constant trajec- 
tory length NAt. 

In addition, we compared simulation results for the full theory at (k,y), with the 
output of a quenched simulation at the corresponding effective coupling value obtained in 
a large-y expansion, 

2 ^ ( 1 
- + ° - 

(cf. Eq. (pip). In table (]2|), we give the mean value of several operators as obtained 
on a 4 3 lattice at k = 0.693, y = 10.0 and in the quenched theory. The agreement is 



feeff = k + 72 + ° ( —4 ) ( 59 ) 
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excellent. Notice that even if the shift in the effective coupling is only 3%, the effects of the 
dynamical fermions can be clearly measured as the observables change quite significantly 
at the critical point k c = 0.693. 

4.3 Phase Diagram 

The phase diagram in fig. [I] was obtained on an 8 3 lattice. As there is no true phase 
transition on a finite lattice, a criterium is needed to locate the phase boundaries. We 
looked for the point where the relevant Binder cumulant equals the value B = 0.8 it has 
at (k = ±0.693 f« k c , y = 0). Since B = 1 deep in the broken phase and B oc 1/L 3 in the 
symmetric one, this provides a clean quantitative criterium which yields a point definitely 
inside the critical region. The width of the critical region decreases as L~ l l u , therefore the 
systematic error in the critical coupling will be at most of order 10 _1 . However, since the 
Binder parameter is a universal quantity, which should stay constant along much of the 
critical lines, the error rather goes as 

L -(u+l/v) ( ie _ 0(io-2)). Thus, this systematic error 
is under control in the full theory as well. We used the standard reweighting method 
to determine the precise location of the points where B = 0.8. 

The total simulation time was 16 days of the 32 Pentium Pro processor parallel com- 
puter RTNN based in Zaragoza. To allow for a correct thermalization, we discarded 100 
integrated autocorrelation times of the relevant susceptibility. This may look utterly con- 
servative, and the MC history indeed seems to stabilize long before that. However, not 
much is known about the exponential autocorrelation time of fermionic algorithms and 
one should be cautious. 

As Eq. (^4|) shows, both at y = oo and at y = we recover the non-linear u-model 
with its well known paramagnetic, ferromagnetic and antiferromagnetic phases. At finite 
y, we expect these phases to extend into the (k,y) plane. In fact one can quite precisely 
anticipate the critical coupling from the strong coupling formula ( |59|) and the quenched 
critical points k^ -00 ^ = ±0.693. Using the reweighting method, the phase transition lines 
can be determined down to y ~ 2.0. In fig. ^ the variation of the Binder cumulant and 
the susceptibility around the two critical couplings is shown for y = 2.0. 

In the small- y region, the effective action up to 0(y 2 ) does not only renormalize k, but 
also introduces additional couplings, due to the non-locality of the matrix K~ l occurring 
in the weak-coupling expansion. Therefore, we do not have an estimate for k as reliable 

Table 2: Comparison of observables in the full theory (]3|) at (k — 0.693, y = 10.0) 
and in the quenched model both at the corresponding value of k eS and at k c = 
0.693. We have 140,000 unquenched trajectories (iV=10, Ar=0.3) on a 4 3 lattice. 
The Metropolis acceptance rate was 65-70%, with an autocorrelation time of 3-4 
trajectories. 



Couplings 




(E) 


d k (E) 


x/v 




fc=0.693 , y= 
k=0.7U , y= 
fc=0.693 , y= 


=10.0 

=0 

=0 


0.4164(6) 

0.41584(14) 

0.3928(3) 


1.134(6) 
1.130(4) 
1.174(4) 


0.3111(7) 
0.3108(2) 
0.2836(4) 


2.378(4) 

2.3779(18) 

2.214(2) 
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Figure 2: Binder cumulant ( |49| ) and non-connected susceptibility ( pf ) as a function 
of k, around the two critical points at y — 2.0. For each critical point, only one 
simulation has been carried out. The other points are obtained with the standard 
reweighting method. 

as in the large- y region fl59|), but we can nevertheless obtain an estimate for k c (y) from 
the MF approximation. We have simulated at several values of the coupling k, for fixed 
y, until the corresponding Binder parameter crossed its critical value. A more accurate 
result for the critical point was later on obtained with the reweighting method. In fig. £3, 
we have plotted the relevant Binder parameter and susceptibility for k values near the two 
critical couplings with y = 0.5. 

In fig. ^| we show the variation of both order parameters and Binder cumulants when 
crossing the FM(S)-EM transition line at y = 1.15. We find a strong change in the stag- 
gered quantities, while the non-staggered ones show a smoother evolution. However, the 
non-staggered order parameter is much smaller than its staggered counterpart. This may 
indicate that, although the non-staggered sector is non-critical {B ~ 1), it will eventually 
undergo a phase transition at lower k. A similar behaviour is found when traversing the 
AFM(W)-EM line at k = —1.6 (see fig. |5|), but now the non-staggered quantities show a 
more pronounced signal. The detailed study of these transition lines (order of the phase 
transitions, critical exponents, etc.) requires a finite-size scaling analysis, which is left for 
future work. This study will be much easier if the transition line is crossed varying k, as 
we lack an analogue of the reweighting method for y. 
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Figure 3: Binder cumulant ( |49| ) and non-connected susceptibility ( |48| ) as a function 
of k, around the two critical points at y = 0.5. The data points are from different 
simulations. 



5 Quasiparticle excitations at the MF level 

In this section we explore the relevant excitations involving fermions, with emphasis on 
the strong-coupling region of our model. 

The small-?/ regime has been studied in relation with the mechanism by which leptons 
and quarks acquire their mass through symmetry breaking in the Electroweak sector of 
the Standard Model. Due to the weak coupling there are no surprises. This situation will 
change dramatically when we consider the strong-coupling region, though. 



5.1 Fermionic excitations in the FM(S) and PMS phases 

At very large y, it is natural to attempt a large-y expansion. This can be achieved after 
carrying out the following change of variables: 

$ = (60) 
= ((/)-T)lp. (61) 

Because of the constraint 4> 2 = 1 and the identity (4> ■ r) 2 = cj> 2 l, this transformation has 
unit Jacobian and its inverse satisfies 

$ = (0 • r) V' • (62) 
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In terms of the new variables (dropping the primes) the action takes the form 

x,fi x,y 

where the fermion kinetic term is the usual lattice kinetic Dirac operator, defined in Eq. 
(|12|). After a further rescaling of the ip fields, the coupling y can be moved to the kinetic 
term, where it appears as 1/y. 

Note that this change of variables ( |6C|j6l| ) was implicitly present in the MF calcula- 
tions of the phase diagram in the strong-coupling region as well (cf. Eqs. ( |30| , j3l| )). This 
transformation is the reason that explains that the model is (partly) analytically tractable. 
The interest of finding reliable analytical approaches to strongly coupled fermion systems 
need not to be stressed. 

The fermion propagator {ipx^y) is given by the expectation value of the inverse fermion 
matrix, which in a large-y expansion becomes 

W*$ v ) = (M-y 1 ) = l-U- -K{ct> ■ t) + - 2 K{4> ■ t)K{4> • r) - . . ) ) . (64) 

\y\y y / xy 

This can be viewed as a sum over paths of increasing length connecting x and y (recall 
that K is a nearest-neighbour matrix). 

In the FM(S) phase, there is a non-zero magnetization v = \ {4>)\- Expectation values 
of products of cj> fields on different sites are replaced by the appropriate powers of v. 
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Figure 5: Binder cumulants and susceptibilities when crossing the AFM(W)-EM 
transition line at k — —1.6. 

Corrections to this approximation as well as contributions from paths visiting a given site 
more than once are of higher order in 1/d and are ignored at the MF level. The series 
( |64D can thus be resummed and one finds a propagator 

1/v 

(^)fm(S) = K + y/v ( 65 ) 

which is that of a fermion with a mass y/v. Note that, since v < 1, this is a huge mass if 
y is large. The propagator for the original fermion, before the change of variables ( |60|j6l| ), 
corresponds to the same physical particle; the only difference is in the wave function 
renormalization. 

In the PM(S) phase, v = 0, so at the MF level the fermion would be infinitely massive, 
or in other words, non-propagating. Beyond this naive MF level, however, a large but 
finite mass will be found. This is due to the next-to-leading contributions to the series 
(|64[). The dominant terms are now those involving the expectation value for the nearest- 
neighbour energy z 2 = {cj> x ■ (f> x+ ^), which is of order l/2d and therefore absent at the MF 
level. The resummation of contributions in (^34|) now leads to a fermion propagator with 
a mass y/z, which is even larger than the mass of the fermion in the FM(S) phase. 

The conclusion of this analysis, which is similar to that in (chiral) Yukawa models 
in the Electroweak theory ]2^| , is that the elementary fermion excitations in the large-y 
region are very heavy (hence essentially non-propagating), and therefore play no role in 
the spectrum of light excitations. This holds even stronger in the PMS phase than in the 
FM(S) phase. 
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5.2 Fermionic excitations in the AFM(S) phase 



Here our point of departure is again the form of the action fl63j), which is tailored for 
studying the large-y behaviour. In the AFM(S) phase, we have a staggered expectation 
value for the <f> field at the MF level, which can be taken in the 3-direction, 

<f>,. :r.,r i ) (GO) 



1 



(with e x = (-l)*i+*a+*s). Hence 



where tpx , i = 1,2 labels the two flavours in tp x . So after the change of variables (|6C| : 61) 
the kinetic operator in (^) is still diagonal in flavour. The only effect of the new variables 
is to change the lattice Dirac operator from (^) to 

ve y T 3 K xy . 

Due to this diagonal structure in flavour space, we can concentrate on one flavour, say 
-i/^ 1 ); the other flavour is obtained by taking —v instead of v. In Fourier space, the kinetic 
term for ip^ is given by 

- ivsfap 6 Pjq ±n, (68) 

where 

s/np = ^Jcr M smp M , (69) 

<>p,q±Tr = <5p M ,(? M +7r mod 2tt • (70) 

So we obtain for the inverse of the MF propagator in the AFM(S) phase, 

Mp, q = -i v s/np <5 Pj(?±7r + y 6 Ptq , (71) 
or, in matrix notation for the subspace of the modes coupled in Eq. ([71]), p and p± (jv, ir, it), 



M P,PH^)-\ iv4np y J' ( 72 ) 

To find the quasiparticle excitations in the AFM(S) phase we diagonalize the fermionic 
part of the action (|72|). One obtains 



S = I ip(p) (y - v s/np) tp(p) , (73) 

where 



p 



i/)(p) = -j= [v (1) (p) + ii> {l) {p + K) 
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or, in position space, 

^ = 7=2 



+ . 

The momentum space propagator corresponding to ((73) is thus 



at \ 1 y + vs/np 

S{p) = =— = -= 5-= — — • 74 

y — v s/np y 2 - v 2 }2\ sm P\ 

Since we are working in imaginary time, one would expect quasiparticle poles in S(p) to 
appear at negative values of p 2 . The unusual relative minus sign in the denominator (|74|) 
therefore does not seem to allow for a quasiparticle interpretation, at first sight. 

However, ( |74|) suggests the possibility of light excitations with a relativistic dispersion 
relation around momenta (±^,±^,±^). To see this, consider the denominator in Eq. 
(fFJ|) for small k^=p^± ir/2: 

y 2 -v 2 Y. sin ' PA = (y 2 - v 2 d) + v 2 £ k 2 + 0(k% (75) 

A A 

where d = 3 is the space-time dimension. As long as we are at large enough y, such that 
y 2 > dv 2 (recall v 2 < 1), this dispersion relation corresponds to a relativistic excitation 
with m 2 = (y 2 — dv 2 )/v 2 , in this naive MF calculation. Several comments are in order: 



1. For v=0, we recover the MF result for the PMS phase: the kinetic term in (|75| ) is 
suppressed. 

2. At the MF level, only for (y 2 — dv 2 ) small enough compared to v 2 these fermionic 
excitations, (t • <t>) if), can propagate easily. Since v 2 < 1, this can only happen for 
y 2 not too large. 

3. These would-be excitations are characteristic of the AFM(S) phase. Let us recall 
that in the PMS phase no light fermionic excitations have been identified at the MF 
level. 



5.3 Light bound states in the PMS phase 

We have seen above that the fermionic excitations in the PMS phase are very heavy. We 
will now show that there are bound states of elementary fermions, however, which are 
light. This is done by means of a MF calculation of the double-chain type [ p4[| . 
Consider the propagator for a fermion pair ip x ipx, 

= ( M *,kr,yXk ~ ( M ^a,;y,X,k M ^yJ ■ (76) 

Here M~ l is the single-fermion propagator, a,/3,\,p are Dirac indices, and i,j,k,l are 
flavour indices. Thus, this propagator is really a 16 x 16 matrix. For the moment we keep 
all these indices as they are; later on we will discuss how pairs of them decompose into 
quantum numbers for the composite state. 
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Figure 6: A typical double-chain diagram, connecting sites x and y. The chains 
are parallel in position space. 



Let us concentrate on the first (M 
of M _1 as before, we find the series 



M~ 



term in Eq. (|7q). Using the 1/y expansion 



(M 



-l 



Ml 



x,/3,j;y,\,k x,a,i;y,p,l 



N,N'=0 



( K 

y V v 



N 



y V y 



N> 



x,a,i; y,p,l I 



(77) 

where we have written 4> as a shorthand for (</> • r). It is clear that only terms with 
N + N' even survive in a paramagnetic phase, due to the cf> — ► —<f> symmetry, thus a 
factor (-1)^+^' has been dropped. Since the matrix K connects nearest-neighbour sites 
only, each term in this series can be seen to represent a product of two paths (chains) of 
lengths N and N' respectively, connecting site x with site y (so, if the "distance" between 
x and y is even(odd), both N and N' will be even(odd)). 

We will attempt to sum the complete series, to leading order in 1/d, where d = 1+2 = 3 
is the Euclidean space-time dimension. For this, we need the spin-spin propagator, which 
in this approximation is extremely short ranged 



1 



■ab 



(78) 



Expectation values of the type (4> x -4> x+ jx) are of order 1/d, and others are suppressed even 
stronger. Thus, assuming (|78|), we observe that any term in the series which contains <f> x 
for a given site x only once or an odd number of times will vanish due to ((f)) = 0. When 
the site is visited twice, it follows from 2 = 1 that the contribution from the fields is 
proportional to \& ab ■ Thus each site along the chains connecting x and y must be visited 
an even number of times. One class of diagrams fulfilling this requirement consists of the 
so-called 'double-chain' diagrams, where the propagation of both fermions between x and 
y follows the same path in position space (see figure As was convincingly argued in 
Ref. [24], this class saturates the dominant diagrams in the 1/d expansion. Indeed, one 
can easily check by concrete examples, how deviations from double-chain behaviour induce 
additional powers of 1/d. We shall also assume that these double chains are self-avoiding 
(this is allowed at first order in 1/d). 
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Figure 7: A matrix-product term contributing to the flavour structure. 



Our task is thus to sum up all double chain diagrams connecting x and y. Let us first 
consider the flavour structure. Using (78) one finds that 



{(<f>x- T )jk(<t> x -T)u) = \^2r°- k T°i = -(-S jk Sii +2SjiS ik ) 



(79) 



^From this and from the ultra- local correlations we are considering (cf. Eq. (|78|)), it follows 
that the product of 2{N + 1) factors of (cf) ■ r) along a double chain of length N visiting 
the points x = xq,xi, ...,y = xjy (cf. Eq. fl77|)) is 



N 



II te 

ra=0 



- 1 x,j;y,k 



N 



n te n , 



n'=0 



(80) 



x,i;y,l' 



To calculate P and Q, it is convenient to represent the general term contributing to the 
above matrix product as in figure |?]. A graph contributing to 5jk5u will have an even 
number of crossings, while diagrams contributing to SjiSik jump an odd number of times. 
Each crossing contributes a factor |, while non-crossings yield factors -§ (cf. Eq. (j7|)). 
Now, P and Q can be easily obtained using binomial summation: 



N 



lite 

,n=0 



x,j;y,k 
1 
3 



N 



n te n , 



.n'=0 



x,i; y,l I 



l\ N 1 



■^(Sjkdii + SjiSik) + (-1) -{5jk5u - SjiSik) , (81) 



where we have separated in a term symmetric under (ji) *-* (ij) and an antisymmetric 
one (this will be needed for separating the contribution to different quantum numbers). It 
is remarkable that the flavour contribution only depends on the double-chain length, but 
not on its shape. This allows for a total factorization between flavour and Dirac indices. 
Next, consider the Dirac structure. One gets products of matrices 



along the double chain, where 



(82) 



(83) 



One readily finds that 



■^xy — ^-^xy-^xy ~ \^y,x+fl + ^y,x—fi)- 



(84) 
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Thus, we need to calculate 



£ 

{/in} 



u 



fin 



(85) 



x,/3,a;y,X,p 



where the sum is extended to all the lattice paths (denoted by {fJ. n }) of- length N connecting 
x and y. Now, we can extend the sum to all length- A lattice paths starting at x, because 
paths not connecting x to y will contribute a zero x y entry. This can be also understood by 
realizing that once the chain has arrived at Xj, there are 2d possible directions to continue 
the chain. These are added up by summing Eq. (]8^) over fi. At the next site, we do the 
same for the next step along the chain. The contributions of all double chains are therefore 
added up when we take the product of these /i-sums along the chain. Corrections due to 
backtracking (2d — > 2d — 1) are down by \/d. 
So we need to calculate powers of the matrix 



(86) 



One way to do that is to write it out explicitly as a 4 x 4 matrix in the space spanned by 
the vectors (/?, A) =(1,1), (2,2), (1,2) and (2,1), in that order. One finds that it equals 



V 



A 3 
A 1 -A 2 





A 1 — A 2 
A 3 





\ 







-A 3 A 1 + A 2 

A 1 + A 2 -A 3 J 



(87) 



It can be diagonalized in this 4x4 space. The eigenvalues, up to the factor 1/4, are found 
to be 



A" 
A 4 



A - 2A" 
-A, 



(^ = 1,2,3), 



where 



.4 



□ + 2d, 



(89) 



(90) 



and □ is the lattice discretization of the d'Alembertian ^Z^d^d^. The A power (see 
flS5"D) of the matrix ( |S6"| ) is now easy to calculate. 

In order to collect the factors and sum up the contributions, let us go back to Eq. 
(f76|). We see that we need to antisymmetrize each term in (Af with respect to the 

simultaneous interchange of a,i with (3,j. This gives a sum of two terms, one symmetric 
in a <-> (5 and antisymmetric in i j (corresponding to a composite state which is a 
Dirac vector and a flavour singlet), and one vice versa (singlet in Dirac space, vector 
in flavour space). Note that Eq. ( |sl| ) has already been written as a sum of symmetric 
and antisymmetric terms. The symmetric and antisymmetric parts of the Dirac structure 
correspond to Eqs. (|88|) and (|89|), respectively. 

Collecting the various factors, we can carry out the geometric sum over in Eq. ( |7"6| ) 
and we find the following propagators for the composite states: 
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• a Dirac vector - flavour singlet with propagator 

-□ + 2A» - 4y 2 - 2d 
where a, b are the Dirac vector indices 

• a Dirac singlet - flavour vector with propagator 

-85 u 
-□ - 12y 2 - 2d 

where /, J are the flavour vector indices. 

These have the form of massive bosonic propagators, up to the following caveat (of 
course, higher order corrections in 1/d may induce shifts in the precise location of the 
poles, as well as their residues). 

The propagators in ([)!]) contain the matrix 2A^ in the denominator. However, this 
term must be ignored since it is sub-dominant in 1/d, compared with the (lattice) d'Alembertian 
□ . 

The numerator of the propagator (^) carries a delta function only, instead of the usual 
tensor structure — d^d u /m 2 . This is also an artifact of the 1/d approximation. 

Notice also that the terms which would play the role of a mass squared in the denom- 
inators have an apparently wrong sign. However, it is easy to check that the composite 
field e x if) x ip x (where e x = (— l)^ 2 ^ as usual) does lead to a massive Dirac singlet - flavour 
vector propagator with mass squared m 2 Q ^ = 12y 2 — 2d= 12y 2 — 6. Similarly, one obtains 

a massive Dirac vector - flavour singlet with a mass squared ^ = 4y 2 — 2d = Ay 2 — 6. 
We thus conclude that the right interpolating field is € x ip x ip x p4|]. 

The conclusion is that we find massive bound states of fermions in the PMS phase. 
They are bound by the strong interactions with the spin waves. These composites are 
lighter than the elementary fermions in this phase, when y moves away from the value oo. 



(91) 



(92) 



6 Conclusions 

In this work we have concerned ourselves with the general features and the analytical and 
numerical study of the lattice model given by expression (Q) . 

From the numerical side, we have developed a new method that exactly solves the 
technical problem related to the length-1 constraint on the spin variable. 

The model describes qualitatively some of the properties of the doped copper oxide 
compounds [|], || and has interesting properties in the strong coupling regime. In fact, 
at the Mean-Field level, no light fermion excitations have been identified in the FM(S) 



and PMS phases. However, in the AFM(S) phase, see section 5.2, light excitations around 
momenta (±f ,±§,±f) have b een found. Its possible relevance for the doped copper 
oxide compounds has been noticed in (l|, g, J|. 

Concerning the PMS phase, (see figure |l|) the situation is also interesting. While the 
fermionic excitations in this phase are very heavy (see section |5.1| ) we have found light 



bound states of fermions (see section 5.3). They are spin singlet bosonic states of charged 
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fermions bound by the strong interactions with the spin waves. A similar result has been 
found in the model of reference [24]. 

The next step is to study the model in the presence of chemical potential and at finite 
temperature (after going to 3+1 dimensions). In fact, as we have proved in section [2|, the 
fermion determinant is still real after the introduction of the chemical potential. 
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